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Abstract 

Protein functions in cells may be activated or modified by the attachment of several 
kinds of chemical groups. While protein phosphorylation, i.e. the attachment of a 
phosphoryl (PO3 ) group, is the most studied form of protein modification, and 
is known to regulate the functions of many proteins, protein behavior can also be 
modified by nitrosylation, acetylation, methylation, etc. A protein can have multiple 
modification sites, and display some form of transition only when enough sites are 
modified. In a previous paper we have modeled the generic equilibrium properties 
of multisite protein modification (R.Chignola, C. Dalla Pellegrina, A. Del Fabbro, 
E.Milotti, Physica A 371, 463 (2006) ) and we have shown that it can account 
both for sharp, robust thresholds and for information transfer between processes 
with widely separated timescales. Here we use the same concepts to expand that 
analysis starting from a dynamical description of multisite modification: we give 
analytical results for the basic dynamics and numerical results in an example where 
the modification chain is cascaded with a Michaelis-Menten step. We modify the 
dynamics and analyze an example with realistic phosphorylation/dephosphorylation 
steps, and give numerical evidence of the independence of the allosteric efi^ect from 
the details of the attachment-detachment processes. We conclude that multisite 
protein modification is dynamically equivalent to the classic allosteric effect. 
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1 Introduction 



Reversible cliemical modifications of proteins are well-known to play a pivotal 
role in the dynamics of the biochemical networks which allow a cell to con- 
vey and translate information from environmental signals to processes such as 
cell activation, proliferation and death pl2] . Moreover, some important bio- 
chemical paths are known to behave as irreversible on-off switches [3j4j, and 
this switch-like character has been associated to the chemical modification 
dynamics of proteins on multiple aminoacid residues or domains [5]. Multisite 
phosphorylation is the foremost example of protein modification because it 
shows up ubiquitously in many important biochemical paths, but in addition 
there are several other multisite modification mechanisms, like acetylation, 
methylation, etc., which act at all levels in the biochemical control networks 
(see [6] for a recent review). Since multisite protein modification (MPM) is 
present in all eukariotes (yeasts, plants, and animal cells), it appears to be an 
evolutionary conserved mechanism that regulates biochemical thresholds and 
switching mechanisms. 

These considerations indicate that MPM is an essential component of many 
biochemical networks. However biochemical networks are complex entities 
where many processes are intertwined with one another and seem to be un- 
approachable with analytical tools: for this reason networks that incorporate 
MPM have been studied numerically in an effort to understand the role of 
MPM itself (see, e.g., [7]). 

Here we attack the problem from a different standpoint: we assume that MPM 
steps have a common character and that they behave much like discrete com- 
ponents in an electronic circuit. Therefore we start by studying MPM in iso- 
lation: in this way it is possible to understand the role of MPM even without 
embedding it in a larger network, and we can produce a few analytical results 
before resorting to numerical methods. 

We have already discussed some of the nontrivial features of MPM in [S] : in 
the following section we briefly review the model and the concepts introduced 
in [8] . Section [3l where we derive the basic set of differential equations, and 
section HJ where we find the equilibrium concentrations, both review previous 
results and introduce the probabilistic interpretation that we have already 
used in [S]. Next we analyze the equations and find the system behavior for 
small deviations from the equilibrium values in sections and O We use the re- 
sults of section [6] to synthesize noise spectra in section [71 We obtain numerical 
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results on the full nonlinear system in section [8] and discuss the general va- 
lidity of the previous results in section [91 We include a downstream catalyzed 
Michaelis-Menten reaction in sections dU] and [TTl We introduce a more real- 
istic phosphorylation/dephosphorylation dynamics in section [T2] and use it to 
produce new numerical results in section [131 Finally we draw our conclusions 
in section [TH 



2 Stochastic analysis of multisite protein modification 

Here we summarize very briefly the results given in [8]: in that paper we 
introduce the reaction scheme shown in figure [H which is very close to that 
considered by Monod, Wyman, and Changeux [H], where a chemical species 
B can modify a number of sites on protein A. We also assume that the N 
modifications sites are all equivalent and that the modification dynamics for 
each site is independent from those of the other sites: this means that we 
consider the states with n modified sites (see, e.g., the transition chain 
shown in figure 2 in [8]). Then, if the single chemical modification dynamics 
is fast with respect to the observation time we can forget the dynamics of 
the transition chain and even the chain itself and concentrate instead on the 
equilibrium probabilities. If the protein becomes activated when the number of 
modified sites is larger than a threshold value rifhr, then from the equilibrium 
probabilities Pn of the modified states A^, we show that the concentration of 
the activated form vs. B has a sigmoid behavior with exponential tails, and 
this defines a very sharp biochemical threshold. The threshold turns out to 
be robust, i.e., it has a reduced sensitivity to parameter changes, which is 
further reduced as and rithr grow. Moreover, when we couple a downstream 
Michaelis-Menten process pUfTT] we find that MPM can produce large delays, 
that once again depend on the number of modification sites N and on rithr, 
and can span several orders of magnitude, thereby providing a link between 
the fast time scale of molecular reactions and the slow pace of cell growth 
and proliferation. We are also able to relate the model to the standard Hill 
phenomenology, which acquires a precise meaning in this context. And yet, 
the approach in [8] is incomplete because a real treatment of the dynamics 
is lacking, and to make further progress we must turn to a better dynamical 
description. 



3 Dynamical model of multisite protein modification 

The model from which we start, and which we modify later to introduce an 
activation threshold, is a classic model in the theory of allosteric activity PITU] . 
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The derivation of the equations is given, e.g., in [TO] section 2.4, and is based 
on the following simplifying assumptions: 



• all sites are equivalent; 

• the occupation of a given site is not influenced by the activity of nearby 
sites; 

• the number of modification sites is constant throughout the process (i.e., as 
the protein modification proceeds and possibly changes the protein shape, 
no new sites are added nor any existing sites are removed); 

• the behavior of the protein depends only on the total number of occupied 
sites, so that the state alone An actually characterizes the protein activity; 

• we assume that the probability of multiple modification events is negligible, 
and therefore we consider only transitions to neighboring states (i.e., there 
are no transitions from An to An+An with |An| > 1). 

• the on-off rates A;+ and fc_ remain fixed and do not depend on modification- 
induced changes. 

From these assumptions one finds the differential system 
d[Ao] 



dt 

d[An] 
dt 



d[A 



N 



dt 



-Nk4Ao][B] + k4Ai] 

-nk^An] - (AT - n)k+[An][B] + {N-n + l)fc+[A„_i] [5] 

+ in + l)k_[An+,] (1) 

-Nk_[AN] + k4Ar,^i][B] 



Summing all the equations we find that A is conserved, and introducing the 
constant value [A]o (total concentration of A, which includes both the unmod- 
ified and the modified forms of A) we write the conservation equation 



E[An] = [A]o (2) 



n=0 



For the moment we also assume that total quantity of B remains fixed, so 
that at any time the following conservation equation must hold as well 

Af 

Y: n[An] + [B] = [B]o (3) 

n=l 

where the square brackets denote the molar concentrations, [B] is the concen- 
tration of the free molecules of B, and [B]q is the total concentration of B 
assuming that all B's are detached. 
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We wish to stress that the nonhnear dynamical system described in this section 
is still highly idealized, it is chemically closed and in thermal equilibrium, and 
is not actually useful until it is coupled with the cellular environment: we 
accomplish this in section [TD], where we relax condition (for a discussion of 
the features of closed and open biochemical systems see, e.g., the introduction 
of the review paper [12j). 



4 Equilibrium values 

As we mentioned above, the differential system ([T]) is in the textbooks, and 
the equilibrium solution is well-known |10|: in this section we review the basic 
results and recast the equilibrium solution in a suitable form. We introduce 
the auxihary variables Pn = [An]/[A]Q, r = {k+/k_)[B] and b = [B]/[A]q, 
and the reduced parameter s = {k+/k_)[A]o, so that r = sb. The p„'s can be 
reinterpreted as the probabilities of finding a protein molecule with n modified 
sites: with this probabilistic interpretation the differential system becomes 
the master equation for the PnS. We also introduce the corresponding barred 
quantities pn, f and b, which denote the equilibrium values, and if we assume 
in addition that the underlying stochastic process (the chain of individual 
modification events) is ergodic, then in the long-time limit these probabilities 
also give the fraction of residence time in each modified state. If the system 
is taken at equilibrium then the derivatives vanish and the differential system 
reduces to a set of algebraic equations. The solution of the system is [lUj: 



Here r is still undefined, and we need yet another equation: we take the con- 
servation equation ([3]), which can be rewritten as 



where bo = [B]q/[A]q, then, substituting the solution (jlj) in ([S]) it is easy to 
show the new conservation condition 




(4) 



N 



^npn + b = bo 



(5) 



n=l 



Nf 



+ b = k 



(6) 



1 +f 



'0 



and then we find a quadratic equation from which we get eventually the equi- 
librium concentration 



(the solution with the minus sign before the square root is unacceptable be- 
cause it gives a negative concentration). The value ([7]) can be translated back 
to the usual notation so that the equilibrium concentration writes 



[BU = \{-{N[A], + {k./k^)-[B],} 

+ ^{N[A], + {k^/K) - [B],Y + A{k./K) [i?]o) (8) 

and the f value that is necessary to evaluate the p's is just f = {k^/k^)[B]f,q. 

At this point it is also interesting to notice that a directly observable quantity, 
the average occupation level, can be associated to the equilibrium values in a 
direct way: 

(n) = ^npn = r-^ (9) 

n=l i^Jo 

A straightforward calculation also yields the variance of the fluctuations close 
to equilibrium 

Nf , , 

varn = (10) 

The equilibrium values [v4„]eq and [-Bjeg depend on the total concentrations 

[A] o, [i?]o, and on the on-off ratio kj^/k_. Whenever {k^/kJ) is large, the 
equilibrium concentration ([8]) can be approximated as follows 

[BU ^ \ - \BM - (iV[A]o - \B\,)\ (11) 

and we see that there is a breakpoint at \B\q = N[A]q, which is the threshold 
of saturation. This behavior is illustrated graphically in figure [2] which shows 

[B] eq VS. [B]q, and where we have set [A]q = lOyuM, which corresponds roughly 
to the concentration of the most abundant proteins in a cell, k^/k^ = 10^ M, 
which is a common value for the on-off ratio (see \10\, p. 56), and we have taken 

= 16 which is the same as the number of putative phosphorylation sites 
for the Rb protein [TSIITU] . Figure [3] shows the corresponding probabilities p„ 
vs. [-B]o, and figure H] shows the relative concentration I]n>nth^[^n]/[^]o where 
rithr is a threshold number of modification sites as in [8] (in this example 
fithr = 10): although the probabilities in figure [3] change considerably when 
the ratio k^/k_ is varied, the relative concentration of the modified A's above 
threshold is remarkably stable with respect to changes of the /c+//c_ ratio. 
Notice that the inverse ratio = 10~^ M is close both to [A]q and [-B]o, 

and for this reason figures [2] and [3] also include the behavior of [B]eq and of 
the pnS for both smaller and larger values of the on-off ratio. We note that 
for large values of k+/k-, i.e. in the case in which the detachment reaction 
is negligible with respect to site modification, the curves are more kinky and 
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change sharply after reaching the saturation threshold. Figure [5] shows the 
modification level (n) and the standard deviation ^/vaFn. Notice that (n) 
grows nearly linearly until the saturation value is reached, at [-B]o = A^[^]o = 
0.16 niM, and that the standard deviation is usually much smaller than the 
average modification level. Similarly, figure [6] shows the average number of 
modified sites in the A^s that are above threshold, i.e., {rit) = J2n>nthr^Pn' 
and the corresponding standard deviation ^var rit of the number of sites above 
threshold: the standard deviation ^var rif is the largest in the vicinity of the 
threshold. 



5 Linearized dynamics 



We have already noted in section [3] that the system is closed, and therefore 
- even though the equations are nonlinear - we can state on very general 
grounds that it must be stable as well p^|IT3] . In the previous section we have 
examined the equilibrium values, but in all possible biological settings it is 
very important to consider the dynamical behavior of the concentrations [An] 
and of the modification level (n) as well, and we turn again to the original 
differential system ([ID, which we rewrite here using the reduced variables: 



k- {-Nrpo +pi} 

k_ {-npn - {N - n)rpn + {N - n + l)rp„_i + {n + l)p„+i} (12) 
k_ {-NpN + rpN-i} 



where we wish to stress that, in addition to the [A„]'s, also [B] and therefore 
also r are time-dependent quantities. If we introduce the deviations from the 
equilibrium values Ap„ = Pn — Pn and recall that b = — Y.n=i '^Pn = b — 
X]^=i^Ap„, we can linearize the system (fT2l) for small deviations: 



dpo 
dt 

dPn 

dt 

dpN 
dt 
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dApo 
dt 

dApn 
dt 



dAp 



'N 



dt 



N 



: k_ I -Ns[bApo -poJ2 ^^Prn] + Api 



m=l 



N 



■ k_ I -nApn - {N - n)s[bApn -PnYl ^^Pn 



m=l 



+ {N -n + l)s[bApn^i - pn-i mAp^] + {n + l)Ap„+i (13) 

m=l J 

: A;_ <^ -NApN + s[bApN-i - Pn-i "^^Pm] \ 

[ m=l J 



In most cases even the linearized dynamics can only be studied numerically 
[T2] . however here we are able to derive the eigenvalues of the system matrix 
[20] , and thus to evaluate all the characteristic time scales of the system. 



6 Eigenvalues 



The generic linearized equation in the differential system f[TS|) is 



1 dApn 
k- dt 



{{N-n + l)fAp„_i -[n+{N- n)f]Apn + (n + l)Ap„+i} 



N 



+s [{N - n)pn -{N -n + l)pn-i] ^^Pn 



(14) 



m=l 



If we temporarily drop the term proportional to the sum Y,m=i i^^Pm, and 
set /c_ = 1, we are left with a system matrix that does not have a definite 
symmetry, but shows a remarkably ordered structure: 



(r) 



-Nf 1 0- 

Nf -1 - (iV - l)f 2 ■ 

{N-l)f -2-(A^-2)f3- 
(iV - 2)f • 



V 



(15) 



In order to find the eigenvalues of A^^ we introduce the matrix U 



{Vn} 



i <k 



^-^^ J > k 



(16) 
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and its inverse 



j <k 



(17) 



Both Uat and its inverse are lower triangular matrices, with a single degenerate 
eigenvalue A = 1. We use Uat to perform a basis change, and we obtain 



{UatAS^^U^i} . ^ = M,, - [{N -k) + {N- k)f] 5,- 



(18) 



which shows that the eigenvalues of A^^-* are A 



[(A^ -k) + {N- k)f] 



{0<k< N). 

The actual system matrix is A = A^^ + A^-*, where A^^ corresponds to the 
dropped term proportional to J2m=i ^^Pm- The elements of the matrix 
are: 

= Nkspo 



0,k 
k 



l<j<N 



k[{N-j)sPj-iN-j + l)sPj.i] 



—ksp 



'N-l 



In the new basis we find 

{u„A<;'ui,'}„^^ 



l<jr<Ar 



{Ua.a!^)u^^} 



+ {N - j)Pj} {N6N,k - SN-i,k) 




and from this result we see that the eigenvalues of the complete system matrix 
are 



A,- 



-k- 




{N-j) + {N-j)sb_ 
l + sb + sE!L-o\N-l)Pi 



(19) 



Vl0<j<Ar-l 
Aat-I 

(where the off rate k_ has been restored to its original value). We have seen 
earlier that the conservation equation ([2]), which translates into the normal- 
ization condition for the probabilities Pn, is automatically satisfied by the dif- 
ferential system (Q, and it can be readily verified that the linearized system 
(fT3l) preserves this condition: this produces the eigenvalue. 

Figure [7] shows plots of the the eigenvalues vs. [-B]o for the example discussed 
in section H] (i.e. with N = 16): apart from the even spacing of the eigenvalues 
Aq to Ai4, one important feature of this plot is the cross-over behavior of the 
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Ai5 eigenvalue. Using the equilibrium concentration b given by equation ([7]) 
and the cross-over condition Xn-2 = ^n-i, after some straightforward algebra 
one finds the cross-over concentration 

_Ns-l 

which approximates to bQ^crossover ~ i.e. [B]q ^ -^[^]o for s 2> 1 (i.e., in 
this case the cross-over value corresponds to the saturation threshold). Notice 
also that for fixed concentrations of A and B, and fixed on-off ratio k+/k-, 
the system moves to higher frequencies (shorter reaction times) as N grows, 
i.e., fluctuations are more effectively damped-off for higher A^'s. 



7 Synthetic noise spectra 



Biochemical reactions where only few molecules are involved are affected by 
molecular noise: this noise often has a deep biochemical meaning [H], and 
is most often studied by Monte Carlo simulation [T^IITS] . However, when the 
eigenvalues of the linearized system are known, as in the present case, it is 
possible to synthesize directly the noise spectra [16]. The eigenvalues deter- 
mine the spectral density of the occupation level of the molecular population: 
because of the random (Poisson) character of the individual molecular events 
the spectral density of the concentrations in the dynamical system, and in par- 
ticular the spectral density S{f) of the modification level n, can be derived 
from the incoherent superposition of the spectral densities of the individual 
relaxation processes associated to each eigenvalue P^|IT7] 

^^^^^ ^ A^ + LfP ^''^ 



We point out that the noise spectrum of the fluctuating modification level has 
a characteristic shape that depends on the concentration of B. The eigenvalue 
distribution for the example discussed at the end of section H] and shown in 
figure [7] indicates that at low concentration [B] the spectral density is roughly 
the superposition of two simple relaxation processes, then close to the thresh- 
old the spectral density collapses to a simple relaxation process (and thus is 
characterized by a 1//^ power-law region). Finally, above threshold the spec- 
tral density has low-frequency white noise region, an intermediate 1// region, 
and a 1//^ high-frequency behavior. The 1/f power-law behavior is limited 
to the range determined by the minimum frequency fmin = fc-( A„)/27r 

0<n<A'' 

and the maximum frequency fmax = k-{ max A„)/27r. 

0<n<N 

All this is further illustrated in figure [S] which shows some synthetic spectra 
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obtained from the eigenvalues shown in figure [71 Figure [8] shows that the 
spectra for the example in section H] have a low-frequency white noise plateau 
and a high-frequency 1//^ noise tail, however when [B]q > N[A]q there is also 
a small 1// noise region which spans approximately one frequency decade just 
as discussed above. However, although the slope change is clear, the 1// region 
is not well defined because of the closeness of the extreme eigenvalues Aq and 

Al5- 



8 Numerical solution of the differential system 

The full differential system ([T]) may be solved numerically with standard in- 
tegration methods: we remark that the (asymptotic) stability properties are 
the same as those of the linearized system, and are guaranteed in our case by 
a theorem due to Poincare and Perron (see, e.g. [20], pp. 161-163), and there- 
fore we do no expect to find any remarkably new features in the numerical 
solutions. We have integrated the differential system ([T]) assuming the values 
at the end of section |H i.e., = 16, A;_ = 1 s~^, k+ = 10® s~^ M~^, and with 
the initial conditions [Ao]t=o = [A]o = lO/iM, [B]t=o = [B]o = 1.2iV[v4]o, and 
[74„]f=o = for n > 0. Figure [9] shows the behavior of the number of occupied 
sites n vs. the dimensionless time variable t ■ k_: the fitting exponential for 
long times is also shown, and the corresponding time constant is in excellent 
agreement with the value computed in the previous section (i.e., the maximum 
nonzero eigenvalue). 



9 Opening up the system 

Up till now we have studied MSM in isolation, and - apart from their util- 
ity in computing the noise spectra - it is natural to wonder if the eigen- 
values can also be useful to understand MSM when the biochemical sys- 
tem is open. Here we consider a straightforward modification, we assume 
that each equation in the differential system ([1]) contains an additional term 
— a„[74„] (we include a minus sign because this additional term is usually dis- 
sipative). With these additional terms the system matrix changes: A^'' — 
A^-* — diagn(ao, . . . , a^), where diagn(ao, . . . , a^) is a diagonal matrix with 
diagonal elements Oq, . . . , oat. Using the transformation matrix Uat and its in- 
verse, it is easy to show that the eigenvalues transform as follows: A„ A„— a„. 
For this reason the previous calculation of the eigenvalues retains its value and 
can be used to estimate the behavior of the system even when it is no longer 
closed and thus when the principle of detailed balance no longer holds. 
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In the next section we turn to a more complex modification of tlie system, one 
wliicli involves the coupling to a downstream reaction. 



10 Switched downstream Michaelis-Menten process 



In this section we consider the following set of coupled reactions 



k+ J n=0,...,N-l 

k2 



n = nthr,--;N 



E + S^ES^E + R. 

ki 

where nthr is the threshold modification level mentioned above, which corre- 
sponds to the onset of release of the secondary enzyme E: the modified species 
An changes to A'^ and releases E. The last reaction is a Michaelis-Menten step 
catalyzed by E which acts on a substrate 5* and produces R. We also make 
the additional simplifying hypothesis that A'^ can no longer take part to the 
modification chain, and can reenter the chain only after reabsorbing E. With 
these assumptions the differential system becomes: 



d[An\ 



dt 



d[Ar^ 



d[Ao 
dt 



n<nthr 



dt 



-Nk^[A,][B] + k_[A{\ 

-nk_[An] - (iV - n)k4An][B] + {N-n+ 1)A;+K_i][fi] 
+ (n + l)A;„K+i] 

-nk^[An] - (iV - n)k+[An][B] + {N - n + 

(22) 



+ (n + l)k.[An+i] - kE+[An] + kEAKWE] 



d[A 



N\ 



dt 



In addition to the equations for the [A.^] 's we must also add the equations for 
the [^^]'s and for the enzyme E: 



dt 



kE+[An]-kE-[A',^][E] 



n>nthr 

m 

dt 



N 



N 



Y: kE^An]- E kE-[A'nm + 



d[E] 



MM 



n=nthr 



n=nthr 



dt 



(23) 
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where the derivative d[E\MM /dt is the contribution of the Michaehs-Menten 
step: 



d[R] 



h[ES] 



(24) 



dt 



d[S] 



h[E][S] + k2[ES\ + s{t) 



(25) 



dt 

d[E]MM 

dt 



h[E][S] + k2[ES] + h[ES] 



(26) 



d[ES] 



k^[E][S]-k2[ES]-h[ES] 



d[E]MM 

dt 



(27) 



dt 



(the term s{t) is the rate with which the substrate S is replenished). As before, 
these equations must be complemented by a conservation equation for [B] 
which now writes: 



We have studied numerically the new modified system: we have taken the 
same conditions as in section [8] and in addition we have set nthr = 10 and 
[A^]i=o = 0. Figure [10] shows the behavior of the number of occupied sites 
n vs. the dimensionless time variable tfc_: the fitting exponential for long 
times is also shown: now the decay constant is much larger than that found 
integration shown in figure [9l i.e., when [B] is above the critical value, the 
dynamical system reacts very quickly to environmental changes. The inclusion 
of the Michaelis-Menten part does not change this behavior, and the approach 
to equilibrium is faster. 



11 Cell-cycle control 

The modified system with the downstream Michaelis-Menten step is remi- 
niscent of the way the cyclin-CDK complex phosphorylates the Rb protein 
which then releases the transcription factor E2F [21], which is an important 
step in the cell-cycle, since it leads to the so-called Gl-S checkpoint [18,.19j. 
Even when we leave aside the enormous complexity of cell-cycle regulation as 
a whole and concentrate on an important detail like Rb protein activity, we 
are still left with a very complicated pattern of biochemical reactions (see, 
e.g., the figure depicting the Rb network in [22]). However the model can be 
further simplified taking only some essential elements from this network (see, 
e.g., figure 8 in [23], which is a good introduction for physicists, see also [2^). 
and in particular we assume that: 



N N 



En[AJ+ n[A'^] + [B] = [B] 
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• the cyclin is destroyed during the cell cycle and is synthesized again after 
mitosis; 

• we neglect the difference between cyclin D and cyclin E; 

• there is plenty of ATP in the environment (cytosol) and we assume that 
the total concentration of phosphoryl groups, i.e., both those in the ATP 
bound to the cyclin-CDK complex and the phosphoryl groups bound to the 
Rb protein, rises roughly linearly as cyclin is produced during the early Gl 
phase p5] (these phosphoryl groups effectively represents B); 

To simulate this process we have integrated numerically the set of equations 
([22])- with the condition 

[B]o = Brt (29) 
where By. is the production rate of available phosphoryl groups and we assume 
that this production rate is small in comparison to the natural relaxation 
rates of the system (i.e., the previously calculated eigenvalues) so that the 
considerations of the previous sections which assume a constant [B]q, i.e.. 
By. = 0, still apply. Notice also that the system is no longer closed and that 
the principle of detailed balance does not hold in this modified situation. 

We keep the conditions that we have already used in the previous numerical 
integrations, and in addition we take: kE+ = 10^ s~^, = 1 M~^, 
fci = 10^ s"^ M-\ k2 = ks = 10^ s-\ the initial conditions [E]t=o = [ES]t=o = 
[5'](=o = [-R]f=o = (these values are in the range of values for Michaelis- 
Menten process commonly found in cells, see, e.g., [11], p. 39). 

Figure [TT] shows [R] for two values of the production rate {Br = N[A]o/T 
and Br = 1.5N[A]o/T, where T = 10^ s): these curves are very similar to 
those that we had obtained in [S] using a simple approximation. Although the 
production rates have a 50% difference the curves are quite close, and this 
suggests that multisite phosphorylation helps making the system robust with 
respect to changes in production rate. 

If we use the point at half maximum as representative of the thresholding 
behavior, we can study the threshold position vs. the synthesis rate Br- This 
is shown by the curve in figure [T21 the curve is well fit by a function with a 
power-law component 

W(x) = a+(6/x)" (30) 
where x = Br/ {N[A]q/T) is the relative production rate, so that the relative 
change of Uhr is 

Atthr ^ a{h/xY ABr 

tthr a+{h/xY Br 

In this numerical integration - which, we wish to stress again, represents a 
realistic and important case - we find a <^ b and a ~ 0.83, therefore 

Atfhr ABr 

— — ^ 0.83—^ (32) 

tthr Br 
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so that there is a shght damping of the fluctuations of production rate, and 
this is an additional factor that contributes to the increased robustness of this 
process with multiple site phosphorylation. 



12 Realistic phosphorylation/dephosphorylation dynamics 



The attachment-detachment dynamics in phosphorylation/dephosphorylation 
chains is actually more complex than the modification dynamics introduced 
in section [3] and used to analyze the example of the previous section. However 
the dynamics can be easily modified to include a more realistic attachment- 
detachment process, like that described in |2l], where phosphorylation and 
dephosphorylation proceed as follows 

protein + ATP ""^^^il^^^ protein-P + ADP 

4- • T~> , u r\ phosphatase „ 

protem-P + — protem + P 

where it is assumed that the reactions proceed in an aqueous environment 
with plenty of ATP. Each reaction is actually a Michaelis-Menten step, and 
if we assume the usual quasi-steady-state approximation [26], we obtain the 
new differential system 



d[Ao] _ ,, fcea,,p[i?(^)][Ao] ^ Ka,AB^^^][A,] 



dt K„,^p + [Ao] Km,D + [^1] 

— —n ■ — — ( A — n) - 



dt K„,^D + [An] Km,p + [An] 



Km,P + [An-l] -f^m,D + [^n+l] 

(33) 



dt K„,^D + [An] ^m,p + [An-i] 



where the kcats and the K^^s are the Michaelis-Menten parameters, the su- 
perscripts P and D denote respectively the phosphorylation and the dephos- 
phorylation step, and [5*^'^^] and [5*^'^^] are the initial (total) concentrations 
of the phosphorylating and of the dephosphorylating enzyme (the kinase and 
the phosphatase in the scheme of ref. [21] )• Notice that now the B's no longer 
depend on the attachment-detachment dynamics, and the differential system 
(1331) seems to be essentially different from the original differential system ([1]). 
However the actual values of the K^s are usually large with respect to the 
expected protein concentrations [21] (for a simple estimate of the protein con- 



15 



centrations inside a cell see, e.g., [8]), and therefore from (l33l) we obtain the 
linear system 



d[Ao] 
dt 

d\M 
dt 



d[A 



dt 



'-Nkp[Ao][B^''^] + kD[A,] 

' -nkn[A^] - {N - n)kp[A^][B^'''^] 
+ {N-n + l)A;p[A„_i][5(^)] + {n + l)kD[An+i] 



(34) 



where ko = kcat,D[B^^^]/ Km,D and kp = kcat,p/Km,p- It is easy to see that 
the system matrix is just like f|T5l) . with the substitution f — >• kp[B^^^]/kD, 
and therefore the eigenvalues are those calculated in section [6], i.e.. 



Xk = -\{N-k) + {N-k) 



kp[B(P^] 



CD 



(35) 



{0 < k < N): thus we see that the more realistic attachment-detachment 
process yields basically the same dynamics in the linear approximation. 



13 Cell-cycle control revisited 



The considerations in the previous section suggest that the inclusion of the 
more realistic phosphorylation/dephosphorylation steps in the cascade of sec- 
tion [TT] should not change appreciably the numerical integration. Unfortu- 
nately, to the best of our knowledge, the actual Michaelis-Menten parameter 
values have never been measured and we have taken simple estimates based 
on related measurements [27f28] . i.e., kcat,p ~ 0.001 s~^, kcat,D ~ 0.0025 s~^, 
Km,p ~ 0.92 /iM, Km,D ~ 0.94 /iM. Using these values, in addition to those 
already listed in section [TT], we have integrated numerically the differential 
system f l33p : even though we have taken K^^s that are not very large with 
respect to the concentrations [An], the results are very similar to those found 
in section [TTl Figure [13] shows a single curve for the product R of the down- 
stream Michaelis-Menten reaction, when we assume that the concentration 
of the Cylin-CDK complex Bp that phosphorylates the pRb protein grows 
linearly in time (this is reasonable, see, e.g., [22]). The B production rate is 
assumed to be 1.39 ■ 10"^^ M s"^: this curve bears a striking similarity with 
those of figure [TTJ Finally [H] shows the threshold time vs. the B production 
rate: again, the curve is very well fit by the function tthr{x) = a + (6/x)°, 
where the exponent is now a = 0.72. The power-law behavior is the same as 
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that found in section [TT] and the exponent is also very close to that found 
earlier; the different value of the exponent is obviously due to the more com- 
plete dynamics. These results indicate that the actual attachment-detachment 
dynamics is not important and that the allosteric effect is independent of its 
details. 



14 Conclusions 



In this paper we have produced a detailed study of the dynamics of multisite 
protein modification, and have analyzed both the equilibrium and the dynam- 
ical properties of the system. This paper follows a previous work [8] where we 
have shown that multisite protein modification may be used by cells both to 
set the time scale of a process - changing it by orders of magnitude - and to 
make it more robust against environmental and endogenous sources of vari- 
ability. Initially we have isolated the attachment/detachment dynamics and 
have explicitly calculated the relaxation rates (eigenvalues) of the linearized 
differential system, and thus the characteristic time scales. The system matrix 
has an interesting, nearly ordered shape and - to the best of our knowledge - 
the eigenvalues that we find here were previously unknown. We have also used 
these rates to compute the synthetic noise spectra, which are relevant when 
the number of molecules is small and discreteness plays an important role. 
We have extended these results with numerical calculations, and we find that 
possible memory effects, that show up as power laws in noise spectra (and 
here we recall that the higher the spectral index, the greater the correlation 
between individual modification events), are suppressed when the modifica- 
tion chain is coupled to a threshold process, and the approach to equilibrium 
is fast. We have considered a process which is very similar to the chain of 
reactions that leads to the crucial restriction point in the cell cycle, and we 
have shown that multisite phosphorylation acts in this threshold sta- 

bilizing factor, that helps reduce individual differences between cells that may 
show up as a different cyclin synthesis rate, and thus stabilizes the duration of 
the cell cycle. These numerical results have been obtained with a grossly sim- 
plified phosphorylation/dephosphorylation dynamics, and for this reason we 
have considered next a more realistic dynamics. We find that the conclusion 
obtained in the section on the linearized dynamics still hold, and moreover an 
explicit numerical integration of the more realistic dynamics yields essentially 
the same results as the simpler bimolecular attachment-detachment dynamics: 
this indicates that the allosteric effect is not actually dependent on the details 
of the modification process. 
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Fig. 1. a. The figure shows schematically the system studied in this paper: the 
modification process is represented by molecules B that react with the sites on A 
with on-off rates A;+, A;_. h. We also assume that when at least n^/i^ sites out of 
the possible sites are occupied, molecule A is activated, and releases an enzyme 
E which catalyses a Michaelis-Menten reaction that converts a substrate S into a 
product R. In this paper we study the dynamics associated to the nonlinear system 
that describes this scheme. In this context we attach a probabilistic meaning to the 
results, we derive noise spectra, and study numerically multisite protein modification 
in conjunction with the downstream Michaelis-Menten step. 
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Fig. 2. Equilibrium concentration [B]f>q vs. \B]q for the example discussed in sec- 
tion m (solid line), for a case with the same parameters but with a higher ratio 
kj^/k^ = 10^'' M (dashed line), and for a case with the same parameters but with 
a lower ratio A:+/fc_ = 10^ M (dotted line). The arrow marks the position of the 
threshold of saturation [B]q = N[A\q. 
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0.05 0.1 0.15 0.2 0.25 0.3 
[BJo (mM) 




[BJo (mM) 

Fig. 3. a). Equilibrium probabilities pn vs. [B]q for the example discussed in section 
m The curves for different p„'s are labeled accordingly, b). Equilibrium probabilities 
for the same parameters but with a higher ratio k+/k- = W^^ M: notice that in this 
case piQ reaches saturation as soon as [B]q reaches the threshold level [B]o = 0.16 
mM. c). Equilibrium probabilities for the same parameters but with a lower ratio 
k+/k^ = 10^ M: in this case saturation is approached much more slowly. 
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0.05 0.1 0.15 0.2 0.25 0.3 
[B]o (mM) 

Fig. 4. Relative concentration Yln>nti [^n]/[^]o of the modified A^s that are above 
a threshold value n^hr {iT-thr = 10 in this example). Solid line: parameter values as in 
the example discussed in section [H dashed line: same parameters but with a higher 
ratio kj^/k^ = 10^'' M; dotted line (lowest curve): same parameters but with a lower 
ratio k+/k^ = 10^ M. 




0.05 0.1 0.15 0.2 0.25 0.3 
[BJo (mM) 

Fig. 5. Solid line: average number of occupied sites (n) vs. [B]q for the example 
discussed in section |H dashed line: standard deviation of the number of modified 
sites ^var n. The thin dotted line shows the saturation value, n = = 16 in this 
case. Notice that the average (n) grows linearly until saturation, while the standard 
deviation is almost always much smaller than the average, and decreases for high 
values of the concentration [B]q . 
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Fig. 6. Statistics of the number of occupied sites nt above a given threshold value 
nthr [nthr = 10 in this example). Solid line: average number of occupied sites (nt) 
vs. [B]q; dashed line: standard deviation ^var nt- The thin dotted line shows the 
saturation value, n = = 16 in this case. 
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Fig. 7. Eigenvalues |A„|/A;_ vs. [B]q for the example discussed in section H] (with 
N = 16): the eigenvalues from Aq up to Aat = A14 are evenly spaced, while the 
eigenvalue A15 is the highest for low [B]q, crosses over the distribution of the other 
eigenvalues, and ends up as the lowest eigenvalue for high [B]q. 
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Fig. 8. These figures display some synthetic spectra of the number of n{t) of oc- 
cupied sites, obtained from the eigenvalues |A„|/A;_ vs. [B]q for the example dis- 
cussed in section H] (with = 16) and shown in figure [71 The amplitude scale is 
in arbitrary units and may change for different spectra, a) spectra obtained with 
[B]o = 0.1iV[A]o, [B]q = N[A]q, and [B]q = 2N[A]q (solid lines) and an ideal l/f 
noise spectrum (dotted line): the first two spectra have a low- frequency white noise 
plateau and a high-frequency 1//^ noise tail. The {B\q = 2N[A\q case also shows 
a limited 1// noise region, between the arrows (which mark the position of the 
extreme eigenvalues Aq and A15). b) Close-up of the 1// noise region for the case 
[B]q = 2N[A\q; here the dotted line is an ideal 1// spectrum and the arrows mark 
the positions of the extreme eigenvalues, c) A larger value [B]q = \{)N[A\q moves 
the 1// noise region to higher frequency; once again the dotted line is an ideal 1// 
spectrum and the arrows mark the positions of the extreme eigenvalues. 
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Fig. 9. Behavior of n{t) from the numerical integration of the differential system ([T]) 
(solid curve) with the conditions specified in section [51 The dashed curve is a fit of 
the tail for t > 0.2s with the function ci — C2exp(— At): we find A = 41.5851/A;_, 
which is very close to theoretical value of the maximum eigenvalue (i.e. largest 
absolute value: |Ai5| = 41.5812//i;_). The number of occupied sites asymptotically 
approaches the equilibrium value, here (n) ~ 15.5709, which is slightly smaller than 
the saturation value = 16 (thin dotted line). 
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Fig. 10. Behavior of n{t) from the numerical integration of the differential system 
in section \T0\ without the inclusion of the Michaelis-Menten reaction (i.e., only the 
equations for the A^s and for the ^^'s are included) and with the parameters 
specified in section [TTJ The dotted curve is the fit of the tail for t > 0.022s with 
the function ci — C2exp(— At): we find A « 550 Hz, which is much larger than the 
decay constant found in the example shown in figure [H The number of occupied 
sites quickly approaches the threshold value rithr = 10- 
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Fig. 11. This figure shows the concentration [R] of the product of the downstream 
Michaelis-Menten reaction of section [TTt when we assume that the concentration 
of the enzyme B that modifies the substrate A grows linearly in time. Curve a 
has been obtained with the parameters given in section [10] and with the B pro- 
duction rate Bj. = N[A\()/T, while curve b has been obtained with the higher rate 
Br = 1.5N[A]q/T. We define a threshold level that is 50% of the saturation level 
of R and we find the corresponding times ti and t2- The large (50%) change in 
production rate leads to a smaller change in threshold time. 
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Fig. 12. Threshold time (defined in figure [TT]) vs. B production rate, from the 
numerical integration of section [TTJ The curve is very well fit by the function 
tthrix) = a + (6/x)", where x = B^ / {N[A\q/T) is the relative production rate, 
with a = 81.7 s, 6 = 4.32 • 10^ s^/", and a = 0.83 (fit and numerical results are so 
close that they are indistinguishable from each other in this figure). 
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Fig. 13. This figure shows the concentration [R] of the product of the downstream 
Michaelis-Menten reaction of section [T3l when we assume that the concentration of 
the Cylin-CDK complex Bp that phosphorylates the pRb protein grows linearly in 
time. The parameters are given in section [T3l and the B production rate is 1.39-10^^^ 
M s-^ 
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Fig. 14. Threshold time (defined in figure [TT]) vs. B production rate, from the 
numerical integration of section [131 here the production rate is x ■ (1.39 • 10^^^ M 
s~^). Again, the curve is very well fit by the function tthr{x) = a + (6/x)", where 
the exponent is a = 0.72 (once again, fit and numerical results are so close that 
they are indistinguishable from each other). 
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